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ABSTRACT 

A numerical and experimental study of the growth of 
succinonitrile (SCN) using a horizontal Bridgman furnace 
and transparent glass ampoule was conducted. Two experi- 
ments were considered: one in which the temperature profile 
was fixed relative to the ampoule (no-growth case); and a sec- 
ond in which the thermal profile was translated at a constant 
rate (steady growth case). Measured temperature profiles on 
the outer surface of the ampoule were used as thermal bound- 
ary conditions for the modelling. The apparent heat capacity 
formulation combined with the variable viscosity method was 
used to model the phase change in SCN. Both 2-D and 3-D 
models were studied and numerical solutions obtained using 
the commercial finite element code, F1DAP 1 . Comparison 
of the numerical results to experimental data showed excel- 
lent agreement. The complex 3-D shallow-cavity flow in 
the melt, differences between 2-D and 3-D models, effects 
of natural convection on the thermal gradient and shape of 
the solid/liquid interface, and the sensitivity of simulations to 
specific assumptions, are also discussed. 

1. NASA does not endorse commercial products. Details 
about the products named in this paper were included for com- 
pleteness and accuracy. No endorsement or criticism of these 
products by NASA should be assumed. 

INTRODUCTION 

The manufacture of electronic materials and optoelec- 
tronic devices demands high-quality crystals. It is known 
(Brown, 1988) that the quality of crystals grown from the 
melt is strongly influenced by the interaction between heat 
and mass transport and fluid flow during the solidification 
process. Consequently, vigorous experimental and numeri- 


cal research activities have been conducted in recent years 
in order to enhance the understanding of transport phenom- 
ena and fluid motion during crystal growth (Sparrow, et al., 
1979, Chang and Brown, 1984, Pimputkar and Ostrach, 1981, 
Langlois, 1984, Glicksman et al., 1986, Yeoh et al., 1990). 

Horizontal Bridgman growth is a widely used technique 
in crystal growth research. During horizontal Bridgman 
growth under terrestrial conditions (1-g), the so-called 
"shallow-cavity" convective flow can be quite strong (Arnold 
et al., 1991, Chait, 1990). Generally, shallow-cavity flow is 
complex, three-dimensional in nature, and about two orders 
of magnitude stronger than convection during vertical Bridg- 
man growth. Previous studies show that this gravity-driven 
convective flow can have a significant impact on growth rate, 
solid/liquid (s/1) interface shape and segregation of impurities 
(Brown, 1988, Sparrow, et al., 1979, Chang and Brown, 1984, 
Pimputkar and Ostrach, 1981, Langlois, 1984, Glicksman et 
al., 1986, Yeoh etal., 1990). 

Among the materials used for understanding solidifica- 
tion and crystal growth behavior, a transparent plastic material 
called succinonitrile (SCN) has been gaining increasing in- 
terest from experimentalists (Mennetrier et al., 1991, Chopra 
and Glicksman, 1988, Inotomi et al., 1993). The main ad- 
vantages of SCN for laboratory study of crystal growth are 
as follows. First, SCN is widely used as an analog to metals. 
For example, a SCN-acetone alloy can be used as a model 
material for metallic solidification study. Second, because 
of its transparency, the fluid motion can be simultaneously 
observed through a transparent ampoule. Thus some other 
important physical characteristics, such as the growth rate, 
interface shape, strength of convective flow and even the con- 
centration profile, can be quantitatively determined using op- 
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tical devices. Third, it has a conveniently low melting point 
(T m = 58.24°C) and grows non-faceted. Finally, its physical 
and chemical properties have been well established (Chopra 
etal., 1988). 

In previous publications (Mennetrier, 1991, Yeoh et al., 
1992, Yao and de Groh, 1993), the authors and their co- 
woikers have reported the effects of natural convection on 
interface shape during horizontal Bridgman growth of SCN. 
Fairly good agreement was achieved among experiments and 
two numerical simulations, one based on the finite element 
method (FEM) and a second used the finite difference method 
(FDM). However, those works were restricted to study of the 
no-growth case. In addition, the study of interface shape 
was limited to a mid-center vertical plane (symmetry plane), 
and no direct temperature measurements were available from 
inside the ampoule. 

In the present work, we consider two new experiments 
(de Groh and Lindstrom, 1994). They are a refined no-growth 
experiment with temperature being measured inside the SCN 


sample and a steady growth solidification experiment with 
non-zero growth rate. In the steady growth case experimen- 
tal conditions were measured after steady state was achieved. 
However, the numerical simulation also models the transient 
between zero and steady growth. For the two experiments, 
the interface shape was measured not only at the mid-center 
vertical plane but also at other two-dimensional planes, which 
together provide a three-dimensional picture of the s/1 inter- 
face. Measured temperature profiles on the outer surface of 
the ampoule are imposed as thermal boundary conditions for 
the FEM modelling. 

The primary objectives of this work are as follows: to 
determine quantitatively how conduction and convection in- 
fluence interface shape, to evaluate the numerical model and 
computer code through comparison with experiments, to ad- 
dress some relevant issues in the numerical modelling, such 
as the difference between 2-D and 3-D models, and to provide 
benchmark numerical and experimental data for researchers 
in this field. 


CELL STAND 



Figure 1 . Schematic of the experimental apparatus, z = 0 is at the center of the gap. 


EXPERIMENT PROCEDURE 

The experiments using SCN were performed in the Low 
Temperature Directional Solidification Furnace (LTDSF) at 
the NASA Lewis Research Center; the experimental proce- 
dure has been described in detail by de Groh and Lindstrom 
( 1994) and Yeoh ( 1 992) . A schematic of the experimental ap- 
paratus is shown in Figure 1. The LTDSF is a Bridgman type 
furnace employing two copper jackets, each of which has an 


associated constant temperature bath. The experiments were 
conducted with the furnace in a horizontal orientation and at 
furnace translation rates of 0 pm/s (for non-growth) and 40 
pm/s (for steady solidification). 

The heating and cooling jackets have a 1.1 cm square 
hole into which the ampoule fits. The borosilicate glass am- 
poules used have an outer square cross section of 0.775 cm 
and are 15 cm long with a wall thickness of approximately 
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0.0925 cm. The comers of the ampoule are slightly rounded as 
depicted in Fig. 2(a). The ampoule was filled under vacuum 
with SCN which was purified to eliminate solutal convection. 
The physical properties of SCN and of the glass ampoule have 
been listed elsewhere (de Groh and Lindstrom, 1994). In this 
paper we shall consider two experiments: the no-growth ex- 
periment (referred to as SCN-4) and the steady solidification 
experiment (SCN-8). (In addition to these two experiments, 
woik by de Groh and Lindstrom (1994) also includes a steady 
melting experiment. Because of the limited space, this third 
experiment will be analyzed in a future publication). 

In all experiments, the heating and cooling jackets were 
maintained at a temperature difference of 63°C' with a hot 
zone temperature of 78 °C for no-growth and 75°C during 
solidification. The temperature distributions on the outer sur- 
face of the ampoule were measured with type K (Chromel- 
Alumel) thermocouples attached to the top, bottom, front and 
comers of the ampoule as shown in Fig. 2(b). The tempera- 
ture of the ampoule was influenced by the furnace set point, 
by convection and conduction in the sample, and by natural 
convection in the ambient air environment around the am- 
poule. These effects caused the temperature distribution to 
vary among the top surface, the vertical sides and the bottom 
of the ampoule. 

To measure the temperature distribution, the furnace was 
moved relative to the ampoule in discrete steps during no 
growth. After each step movement of the furnace, the sys- 
tem was allowed to equilibrate for about 10 minutes before 
temperature measurements were taken. This gave the tem- 
perature distribution of the ampoule as a function of distance 
from the s/1 interface. The temperature and position measure- 


ment accuracies are estimated to be approximately ± 1 ° C and 
±0.5 mm. 

In the no-growth case, temperature was also measured 
inside the SCN sample with one thermocouple located at 
the center and the other between the center and die upper 
wall. These two thermocouples, designated in Fig. 2(b) by 
T e and T u respectively, provided temperature distributions 
along two longitudinal directions, i.e. the lines (0, 0, z) and 
(1.54mm, 0, z). 

A microscope and camera were used to examine the s/1 
interface. The interface location and shape was quantitatively 
analyzed at the mid-center vertical plane (MCP), upper-center 
horizontal plane (UCP), top surface (TS) and the front surface 
of the SCN sample (FS). In the coordinate system shown in 
Fig. 2, the MCP corresponds to the (a;, 0, z) plane. UCP is 
in the yz plane at x = 1.5mm for SCN-4 and x = 1.0mm 
for SCN-8. TS represents the intersection between the S/L 
interface and the inner top wall, namely the interface shape 
in the plane. FS refers to the intersection between the s/1 
interface and the inner front ampoule wall, the interface at 
the (2.95mm, y, z) plane. Photographs were taken of the 
s/1 interface only when steady state conditions were reached. 
The photographs were then traced on an electromagnetic dig- 
itizer. The shape of the interface was obtained from the dig- 
itized data. To minimize the inconsistency of the measured 
temperature and interface positions caused by experimental 
uncertainties, a simple averaging technique was used to adjust 
the raw experimental data. These adjusted temperature and 
interface position data are listed in de Groh and Lindstrom 
(1994). 
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(a) (b) 

Figure 2. Schematic of the cross-section of glass ampoule and arrangements of thermocouples, (a) The 
rounded-corner ampoule actually used in the experiment and definition of co-ordinate system, z = 0 is defined 
as the center of the gap between hot and cold zones and the positive z-axis is directed towards the solid. 
6 = 0.0925cm is the thickness of ampoule wall, (b) The simplified sharp-corner model used in the simulation. 
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The Prandtl number, Pr = v/a, for SCN is approx- 
imately 23, where v = 2.6 x 10 ~ 6 m 2 /s is the kinematic 
viscosity and a is thermal diffusivity. The thermal Grashof 
number is given by: 

v * 

where fix — -8.1 x 10 -4 /A, is the coefficient of ther- 
mal expansion. Since we are considering horizontal growth, 
AT = ISA'' is the temperature difference between the hot 
zone and the melting temperature, and the distance over 
which this gradient acts is l c , the characteristic length of the 
system (l c — 2 cm). Under these conditions at unit gravity, 
g = 9.8m/s 2 , the thermal Grashofnumberis about 1.7 x 10 5 . 

NUMERICAL MODELING 
Governing Equations 

The mathematical model used in the present work con- 
siders heat transport in both SCN and the glass ampoule, fluid 
motion in the melt, and phase change at the S/L interface. 
The liquid SCN is assumed to be Newtonian and its motion 
is described by the following Navier-Stokes equation: 

Po + u- Vuj = - Vp + V ■ [//(Vu + (Vu) T )] 

+ />og[l - Pt(T - To)] (1) 

where u is the velocity vector, t is time, po is fluid density, 
p is pressure, p. is viscosity, T is temperature, T 0 is a ref- 
erence temperature, g is the acceleration of gravity, (1 is the 
volumetric expansion coefficient and the Boussinesq model is 
adopted to approximate the buoyancy force caused by density 
variation with temperature. The incompressibility condition 
for liquid SCN is given by 

V • u = 0 . (2) 

The heat transport is controlled by the balance of thermal 
energy 

P° Cp + U ' VT ) = V ' («VT) (3) 

where c p is specific heat, k is thermal conductivity. To model 
the borosilicate glass ampoule we treat the ampoule walls 
as solid and consider only the heat conduction equation (3) 
within the wall regions. Since the temperature levels consid- 
ered were low, radiative heat transfer was neglected in this 
work. 

The boundary conditions associated with eqn. (1) are 
the no-slip conditions imposed at the SCN/ampoule interface. 
For eqn. (3), the interpolated temperature profiles based on 


the experimental measurements are specified on the outer 
surface of the ampoule. Note that the numerical simulation 
incorporates conduction through the ampoule. 

At the s/1 interface, the following phase change condi- 
tions need to be satisfied: 

T 3 [S(t),t} = T l [S(t),t] = T m , (4) 

(k.VT. - k/VT z ) • n = pL~ . (5) 

The subscripts, s and /, refer to the solid and liquid regions, 
respectively; S(t) is the spatial position of the S/L interface; 
T m is the melting temperature; n is the unit outward normal to 
the interface (pointing from liquid to solid); and L is the latent 
heat of fusion. Since at the macroscopic level the interfaces 
studied are flat, we shall not consider the influence of inter- 
face curvature on melting point (Flemings, 1974) and assume 
solidification to occur at the equilibrium melting temperature. 

The FEM Model 

In the horizontal Bridgman configuration, the longitudi- 
nal axis of the furnace is perpendicular to the gravity vector 
as shown in Fig. 3. Since it is reasonable to assume that the 
heat and flow fields are symmetric with respect to the middle 
center plane, only half of the ampoule is modelled in the 3-D 
Simulation. We consider two ampoule cross-sections. One 
has rounded comers, which is closer to the real ampoule used 
in the experiment, while the other has sharp comers. 

In selecting the length of the computational domain, we 
considered only the central 6.25 cm of the total 15 cm length 
of the ampoule. Beyond the region modeled, the experimental 
data indicated that the temperature variations are very small 
and should not affect the solution near the interface. 

The 2-D FEM mesh was built with the 4-node bilinear 
element, in which the velocity and temperature were approx- 
imated by bilinear shape functions and the pressure was ap- 
proximated as piecewise constant The typical 2-D mesh used 
in the modelling has 954 bilinear elements and 1026 nodes. 

For the 3-D FEM model, we used the 8-node linear brick 
element, in which the velocity and temperature are assumed 
to be trilinear and the pressure to be piecewise constant. The 
FEM meshes were generated by FIMESH, a mesh generator 
provided in FTDAP. A total of four 3-D meshes were cre- 
ated for our computation. For each of the sharp comer and 
rounded comer models, we had a corresponding fine mesh 
and a coarser mesh. The two fine meshes are shown in Fig. 4. 

In this work numerical solutions were obtained using 
the FEM program FIDAP, a general fluid dynamics analysis 
package (Engleman, 1993). For steady-state problems, the 
discretized FEM equations can be reduced to a set of nonlinear 
algebraic equations in the following matrix form 

K(U)U = F (6) 
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Borosilicate Glass Ampoule 

I 


Liquid SCN 


interface 


Figure 3. Schematic diagram of the simplified 2-D FEM model for the directional SCN 
growth experiment. For the 3-D model this diagram represents the mid-center vertical plane, 
i.e. the xz plane at y = 0. 



Figure 4. The two finer3-D meshes used. On the left is the rounded corner mesh. This mesh 
has 1 8090 8-node linear brick elements with a total of 20468 nodal points. On the right is the sharp 
comer mesh built with 19872 elements and 22750 nodes. 
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where K is the global system matrix, U = (u ,p,T) is the 
global vector of unknowns (velocities, pressure and tempera- 
ture), and F is a vector that includes the effects of body forces 
and gradient type boundary conditions. For the problem we 
considered, there is strong coupling between the momentum 
and the energy equations through the buoyancy term as well 
as through the convective terms. 

The Enthalpy Method 

To model the phase change in SCN, we used the en- 
thalpy method (Crank, 1984). As a fixed domain approach, 
the enthalpy method avoids tracking the s/1 interface by incor- 
porating the Stefan condition (5) into the following definition 
of enthalpy: 

H{T)= f T pc p (r)dr + pf(T)L (7) 


where T re f is an arbitrary reference temperature and / is 
the local liquid volume fraction (Crank, 1984, Voller and 
Swaminathan, 1990). For an isothermal phase change, / is 
given by the Heaviside step function 



T<T m 

T>T m 


( 8 ) 


where k represents iteration number. 

One of the difficulties inherent to the AHCM is the selec- 
tion of parameter e involved in the calculation of Cf or H e . 
Detailed studies about the effect of e value on solution accu- 
racy and convergence with mesh refinement can be found in 
Shamsundar (1978), Yao and Chait (1993a), (1993b). For the 
AHCM with the linear approximation given in (9), numerical 
tests in Yao and Chait (1993a) shows that the accuracy of so- 
lutions varies with the value of e and there exists an optimum 
e at which the best accuracy may be obtained for a particular 
spatial and time discretization. Unfortunately, the analytical 
prediction of the optimum £ is not available in practice. The 
general rule recommended by Bonacina (1973) is to select 
the value of £ so that the transition zone embraces one or two 
complete elements. This rule was used for selecting e values 
in our computations. Our study suggests that this rule pro- 
duces results that are very close to the optimum e. In enthalpy 
method, the s/I interface is usually recovered by computing 
its isotherm; its resolution is determined by the value of £. 
For the no-growth case the resolution of interface is about 
2e = 0.1°C which is smaller than the estimated accuracy of 
measured temperature. 


In numerical solutions, H as given in (7) usually needs 
to be smoothed. The following linear approximation (Yao 
and Chait, 1993) is used in our computation: 


He(T) = { 


( PsCpT , T < T ro - e 

PsCp(T m — e) + T m — £<T<T m +e 
PsCp + Pic l p Pe L 


+ 


+ ^)[T-(T m -e)}, 


PsCpT m + piC (T T m ) + PaT 


\ T >T m + e . 

( 9 ) 

In this approximation, the latent heat is released over a small 
temperature interval [T m -e,T m +e}, and the value of e cor- 
responds to the half-length of a transition zone. 

Based on the definition of enthalpy an apparent heat 
capacity (AHC) can then be defined as (Crank, 1984) 


fiTT 

C a (T) = — = pc p + pL6(T - T m ) (10) 


where 6 is the Dirac function. Numerical schemes based 
on (10) are often referred to as the apparent heat capacity 
method (AHCM). We use the following temporal averaging 
scheme (Morgan et al., 1978, Dantzig, 1989) to approximate 
the AHC: 




H e (T k ) - HejT*- 1 ) 

j<k _ yfc-i 


( 11 ) 


At the s/1 interface or in the transition region, the fluid 
velocity diminishes to zero. Clearly then, the enthalpy for- 
mulation must be able to account for this velocity behavior 
in the vicinity of the phase-change front. There are a number 
of methods available in the literature (Voller, 1987). In our 
computation, a so-called variable viscosity method (Voller et 
al., 1987,Gartling, 1980) is used to model the velocity behav- 
ior. In this method, the liquid SCN is treated as a Newtonian 
fluid whose viscosity is a function of temperature and takes 
an artificially high value when T is below the melting point. 
This treatment has the effect of immobilizing the solid portion 
of the SCN so that any predicted velocities in the solid region 
are negligible. 

Implementation of the variable viscosity method is very 
simple. It requires only the specification of a viscosity- 
temperature curve. However the artificially introduced large 
discontinuity in the viscosity may cause some numerical dif- 
ficulties in convergence of the nonlinear iteration, especially 
when the segregated solution approach is used and the convec- 
tion is strong (as in this work). Consequently, an incremental 
solution procedure is used. In this solution process, we be- 
gin with a viscosity jump across the interface of about 10 4 . 
Then we restart from the previous solution each time with a 
viscosity jump increase of 10 or 10 2 . Our computation stops 
when a viscosity discontinuity of 10 8 is reached, since our 
experience shows that the solution accuracy is good enough 
at this point. 
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Segregated Solution Approach 

In general there are two distinct solution approaches 
available for solving the nonlinear system (6). The first 
approach solves the whole system in a simultaneous cou- 
pled manner, while the second approach solves each equation 
separately in a sequential segregated manner. Fully cou- 
pled solution approaches, such as successive substitution and 
Newton-Raphson, usually require formation of the global sys- 
tem matrix which includes all the unknown degrees of free- 
dom. While this strategy is cost-effective for most 2-D prob- 
lems, the peripheral storage required for 3-D problems can 
become excessive. 

In contrast to the fully-coupled approach, the segregated 
solution algorithm (Haroutunian et al., 1991a and 1991b) 
avoids the direct formation of a global system matrix. Instead, 
the global matrix is decomposed into smaller sub-matrices, 
each governing the nodal unknowns associated with only one 
conservation equation. These smaller sub-matrices are then 
solved in a sequential manner. As the storage of the individual 
sub-matrices is considerably less than that needed to store 
the global system matrix, the storage requirements of the 
segregated approach are substantially less than that of the 
fully coupled approach. This makes it possible to solve large- 
scale 3-D simulation problems on a workstation. All the 
3-D computations presented in this work were done on a 
SPARCstationf. 

f SPARC is a registered trademark of SPARC interna- 
tional, Inc. 

RESULTS AND DISCUSSION 

In this section we present our numerical simulation re- 
sults and compare them with experiments. 

No-Growth Experiment (SCN-4) 

2-D FEM Modeling 

It is natural to start from a simplified 2-D model which 
considers the mid-center vertical plane of the ampoule (MCP) 
only. The measured temperature profiles given by the ther- 
mocouples T t and T& are imposed on the top and bottom 
boundaries of the 2-D model. Since the resulting global sys- 
tem (6) is much smaller than that in the 3-D case, the fully 
coupled solution approach is used for solving the nonlinear 
system. 

The numerical results suggest that the heat and flow 
fields predicted by the simplified 2-D model are fairly good 
for the no-growth experiment. For example, the 2-D solution 
of temperature along the center z-axis, T c , is very close to 
the 3-D solution and to experimental data, as shown in Fig. 5. 
The S/L interface shape of the 2-D model also compares well 
with the 3-D solution and the measured data, as can be seen in 
Fig. 6. Ourresults indicate that solutions of the 2-D model can 
provide reasonable approximations for temperature, velocity 
and interface shape on the symmetric planes (MCP) of the 
ampoule. 


Visualization of 3-D Flow In the Melt 

The information obtained from the 2-D model is limited 
to the MCP. Although the 3-D model requires much more 
computational effort and resources, it can provide much more 
detail on the interface shapes, transport phenomena and fluid 
flow in the whole domain. The flow pattern on the MCP indi- 
cates that the dominant flow is the so-called shallow-cavity 
Sow (Arnold et al., 1991), which agrees with the experimen- 
tal observation and literature. The primary characteristic of 
shallow-cavity flow is the single recirculating cell in the ver- 
tical (xz) planes as shown in Mennetrier et al. (1991), Yeoh 
(1992), and Yao and de Groh (1993). This strong convective 
flow cell forces hot liquid to flow along the top wall, raising 
and homogenizing the top wall temperature as well as direct- 
ing hot liquid against the upper portion of the interface. The 
above mentioned fluid motion can be visualized using a trace 
of a particle path as shown in Fig. 7, in which a particle is re- 
leased at the point x — 0.26, y = -0.1 &z = -0.6. As seen 
in the side view in Fig. 7(a), the liquid particle turns downward 
when confronted by the interface. After being cooled by the 
interface and the imposed thermal gradient, it returns along 
the bottom wall of the ampoule. The isometric perspective 
view in Fig. 7(b) illustrates the complex three-dimensional 
particle trajectory in the space between the MCP and the rear 
wall. Note that the spatial position of the particle in Fig. 7 
is plotted at each time step. Since the time step is fixed, the 
distance between two particle symbols actually indicates how 
far the particle travels in one time step. Therefore it is clear, 
by examining the distance between the circular symbols on 
the path, that the particle flows much faster near the interface 
than at the other end of the domain. 

Another important feature of the 3-D flow field is sec- 
ondary flow in the planes perpendicular to the z-axes (i.e. the 
xy planes). A visualization of the secondary flows is shown 
in Fig. 8 in which we plot velocity vectors on different xy 
planes. The plot at z = 0 cuts the interface in the middle and 
captures the liquid SCN on the upper portion only. The blank 
lower portion in the z = 0 plot is the solid SCN. Moving away 
from the interface along the negative z-direction, a flow cell 
is first developed on the upper portion around z = -0.2cm. 
This flow cell grows stronger on the xy planes further away 
from the interface. After reaching the maximum strength, 
the upper cell starts decaying. Meanwhile, a second flow 
cell is formed in the lower part as shown by the z = —0.4 
plot This lower cell rotates in an opposite direction to the 
upper cell and is stronger in planes deeper in the hot zone. 
At around z = -0.6cm (about one seventh of the total liquid 
SCN length) the two counter rotating cells are about of equal 
strength. Then the upper cell becomes weaker and weaker, 
being pushed to the upper comer at z = -3.2 and vanishing 
at about z = -3.9. Atz = -4. lcm, the flow is dominated by 
the primary flow pattern again with liquid coming along the 
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Figure 5. Comparison of 2-D and 3-D FEM solutions of temperature variation of T c along the 
center z-axis, i.e. line (x = 0, y = 0,z), with the experimental data. 
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Figure 6. Simulated s/I interface shapes on the MCP based on 2-D and 3-D FEM models 
and comparison with the no-growth experiment (SCN-4) data. 
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Figure 7. The particle path plotted from time 0 to 220 sec. with a fixed time step At = 1.2 
sec for SCN-4. The particle, labelled as P, is released near the interface and MCP. (a) The side 
view; (b) The isometric perspective view. Note that only the liquid SON part of the computational 
domain is displayed here. 



Z = Ocm Z = -0.2cm J_ Z = -0.3cm Z = -0.4cm Z = -0.5cm 



Z = -0.6cm j_ t Z = -2.2cm Z = -3.2cm Z = -3.9cm . z = ~ 4 - lcm 


Figure 8. Visualization of the secondary flows on different xy planes for the no-growth SCN-4 experiment. 
The plots are based on the finer mesh for the simplified sharp corner model. Note that the plots are viewed 
from negative ^-direction and the hot end of ampoule is z = -4.2cm. 
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Figure 9. Computed no-growth 3-D FEM solutions of temperature variation of T u along 
the upper z-axis, i.e. line (a; = 1.54mm,y = 0 ,z) using the finer round-corner and the 
sharp-corner meshes. 




Figure 10. Modelled S/L interface shapes during no-growth solidification using the round-corner and 
the sharp-corner FEM meshes and comparison with the experimental measurements, (a) On the front 
surface of the SCN sample (FS, y = -2.95mm); (b) The upper center plane (UCP, x = 1.54mm). 
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bottom wall, changing direction and flowing upwards along 
the hot end of the ampoule, then changing direction again and 
flowing back along the top wall toward the interface. 

Effects of Natural Convection 

Our previous studies and the present work indicate that 
the convective flow in the melt has a great impact on the shape 
of s/1 interface. This is shown in Fig. 6 and 12 where the 
interface obtained from the conduction-only solution is flat; 
the interface is highly curved when convection is included. 

Another interesting phenomenon related to natural con- 
vection is the temperature distribution in the neighborhood 
of the interface. As seen in Fig. 5, the temperature variation 
is smooth near the interface (around z — 0) with a relatively 
wide diffusion layer when no convection is considered. How- 
ever, when convection is added, there is a sharp change in the 
slope of the temperature curve and the diffusion layer becomes 
much narrower near the interface. Note the good agreement 
between the simulation with convection and the experimental 
data in Fig. 5. This narrow diffusion boundary layer and the 
dramatic change in slope render a great challenge to numer- 
ical modeling. We believe this phenomenon may be one of 
the main causes of difficulty in attaining convergence. 

As noted earlier, for the cases considered, the thermal 
Grashof number was about 1.7 x 10 5 and the Prandtl num- 
ber is about 23. At first, this Grx seems moderate, at least 
compared to some other horizontal systems examined by de 
Groh and Nelson (1994). However, SCN’s relatively large 
Prandtl number enables larger transverse thermal gradients 
to develop, thus making possible the highly curved interface. 
Materials having a lower Prandtl number tend to result in 
flatter isotherms and interfaces (Gadonniex et al., 1994 and 
Rouzaud et al., 1993). This is reflected in the Rayleigh num- 
ber, Rai - = Gtt x Pr, which was about 4 x 10 6 in this 
study. For two other cases of horizontal plain front solidifi- 
cation using Bi-Sn alloys, which resulted in flatter interfaces, 
Rayleigh numbers were in the range of 4 x 10 5 (de Groh and 
Nelson, 1994). Thus for SCN the buoyant force/viscous force 
ratio is comparatively larger as compared to metals, making 
SCN generally more sensitive to gravity driven convection 
than metals. 

Effects of Ampoule Corner Shape 

To study the effects of ampoule shape on the thermal and 
flow fields, and especially the effects on the interface shapes, 
we examined both round-comer and simplified sharp-comer 
meshes as depicted in Fig. 1 . The comparison between the two 
sets of solutions indicates that the shape of the ampoule comer 
has very little effect on the numerical solution. For example, 
on the MCP, the modelled temperature, velocity and interface 
shape based on the two meshes are almost identical. As a 
typical check, the temperature distribution of T u on the MCP 
in Fig. 10 shows less than 1% difference between the solutions 
of T u from the round-comer and sharp-comer meshes. 

Slight differences are observed only in the vicinity of the 


comers. In Fig. 10(a) we compare the interface shape at the 
front surface of the SCN sample (FS). As shown, the interface 
shape given by the round-comer mesh is slightly closer to the 
experiment than the shaip-comer mesh solution. Another 
comparison of interface shape on the UCP is presented in 
Fig. 10(b), in which only a very small difference can be 
observed near the ampoule wall (at y = ±2. 95mm). 

Steady Solidification Experiment (SCN-8) 

2- D FEM Modeling 

As in the no-growth case, the 2-D model is adopted to 
simulate the middle center plane (MCP) of the ampoule. The 
measured temperature profiles of the thermocouples T t and 
Tb (shown in Fig. 11) are used to specify thermal boundary 
conditions on the top and bottom boundary of MCP. A linear 
interpolation is used when the nodes fall between the raw 
experimental data points. To model the movement of the 
furnace (and hence the transient thermal environment), we 
assume the temperature profiles of T t and T& do not change 
with time, but translate towards the hot end with the same 
constant speed of 0.04mm/sec as the furnace. Since there is 
no direct experimental data available for an initial solution 
we use the steady (no-growth) analysis results as the initial 
condition. In our computation the Euler backward scheme 
with a fixed time step At = 2 sec is used for the time inte- 
gration. The resulting global FEM system is solved by the 
Newton-Raphson iteration scheme. 

To check the validity of the 2-D model we compare the 
solutions of interface location on MCP with experiment and 
present the results in Fig. 12. The conduction-only solu- 
tion does not agree well with the experimental results and 
shows once again the significance of natural convection ef- 
fects on the interface shape. In section 4.1 we have shown 
that the 2-D model provides a very good approximation for 
the interface shape on MCP in simulating file no-growth ex- 
periment (SCN-4). However this conclusion does not hold 
for modelling the solidification with non-zero growth rate. 
Fig. 12 shows the maximum difference between the 2-D so- 
lution (at t = 400sec) and the experiments to be about 3mm. 
This rather large difference is considered unacceptable and is 
much greater than the uncertainty in interface location in the 
experiments, which was estimated to be ±0.5mm. 

3- D FEM Simulation 

The main reason for the inadequacy of the 2-D model 
during growth at the MCP is that the 2-D model does not 
have any knowledge of the 3-D effect from the rear and front 
ampoule walls. This does not seem to be a problem when the 
interface is not moving, such as the case in SCN-4 because 
there is less curvature in the yz plane. However when the in- 
terface moves, the 3-D effects become much more important. 
Our study suggests that a full 3-D model has to be used for 
simulating the growth of SCN-8. 

Because the segregated solution approach is currently 
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Figure 11. The measured temperature distribution of T t and T& along the ampoule axis. The 
imposed T t and T& are obtained by linear interpolation between experimental data points and used 
in the 2-D modeling as the top and bottom temperature boundary conditions. During growth the two 
curves are translated with a steady speed of 0.04mm/sec towards hot end. 



Z(mm) 

Figure 12. Simulated S/L interface shapes on the middle center plane (MCP) based on the 2-D FEM 
model and comparison with the experiment. The conduction-only solution is based on steady analysis with 
no convection. The initial solution is obtained by steady analysis considering both conduction and convection 
but no growth. The transient solution considers conduction, convection and a steady translation of the furnace 
(steady growth). 
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Figure 13. Development of S/L interface shapes on the MCP given by the transient 3-D 
FEM solidification analysis and comparison with experiment during horizontal solidification. 



Figure 14. Simulated S/L interface shapes during solidification on the upper center plane 
(UCP, x « 1.0mm) and comparison with experiment. 
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Figure 15. Computed interface shapes on the top surface of the SCN sample (TS, 
x = 2.95mm) based on 3-D transient analysis for SCN-8 and comparison with experiment. 



Figure 16. Time history of the interface position at the center of top surface (a; = 2.95, y = 
0) and the center (x = 0, y = 0) of SCN sample. The speed of the moving interface is given 
by the slope of the curves. 


14 




9 


» 


* • 


not available for transient analysis in FEDAP, the 3-D results 
presented in this section were obtained using the fully coupled 
solution approach and based on the coarse 3-D mesh with the 
simplified sharp comers. 

The numerical evolution of the interface shape at the 
MCP is shownin Fig. 13. Our results indicate that the solution 
becomes steady after about 200 seconds. The interface shapes 
at t = 200 and 400 sec are quite close to the measured 
data (Fig. 13). The simulated interface shapes on the UCP 
and TS planes are presented in Fig. 14 and 15, respectively. 
The agreement between the 3-D numerical solution and the 
experiment is excellent. 

The time history of two interface positions are shown in 
Fig. 16. One is at the top surface of SCN sample, i.e. along the 
line (x = 2.95mm, y = 0, z) and shown by the curve referred 
to as “surface” in Fig. 16. The other is at the center of the 
SCN sample. The slope of the two S(t) curves represents 
the simulated growth rate, and the vertical distance between 
these two curves indicates the interface deflection between the 
top surface and the center. As we can see from Fig. 16, the 
interface shape is fully developed and the movement of the 
interface becomes steady when t > 200sec. The predicted 
growth rate, which is 0.04 mm/sec, is equal to the furnace 
translation speed used in the experiment. 

CONCLUSIONS 

Perhaps the most obvious advantage of using SCN for 
crystal growth research is its transparency, which allows quan- 
titative determination of the s/1 interface shape and simulta- 
neous observation of important transport phenomena. The 
numerical and experimental investigation presented in this 
paper provides two complete benchmark tests for the hori- 
zontal Bridgman growth of SCN under terrestrial conditions. 
The excellent agreement between the numerical results and 
experimental measurements not only serves as a validation 
of the numerical model but also demonstrates the importance 
of the collaboration between numerical modeling and exper- 
iment. 

For the two cases examined, shallow-cavity convective 
flow is confirmed, by our numerical and experimental obser- 
vation, to be the primary flow pattern. This convective flow 
and its dominance in heat transfer have a significant impact on 
the shape of the s/1 interface. As viewed from the liquid side, 
the s/1 interface in both cases is concave in the upper half of 
the ampoule. This is due to the above-mentioned converting 
flow which brings relatively warm liquid from the hot zone 
and rams it against the upper part of the interface. This flow 
of warm liquid moving toward the interface along the top am- 
poule wall with cooler liquid returning along the bottom wall, 
results in a dominating thermal gradient in the x-direction 
(with the top hotter) and interface asymmetry as viewed in 
the vertical (xz) planes. During solidification, as the longitu- 
dinal thermal gradient is translated along the ampoule in the 


negative x-direction, the liquid SCN gets cooled immediately 
near the ampoule wall and the s/1 interface in contact with the 
ampoule moves with the same translation speed. However, 
the inner regions of the sample must conduct the heat through 
the poorly conducting SCN, causing solidification near the 
ampoule center to lag behind the outer edges. This results in 
a much more deflected (curved) interface for the solidifica- 
tion case. In addition to the effect on interface shape, natural 
convection in the melt has also a significant effect on the tem- 
perature distribution near the interface. Our experimental and 
numerical results indicate that the convection causes a very 
sharp thermal gradient in the liquid near the interface. 

Another interesting aspect of the complex 3-D flow mo- 
tion in the melt is the occurence of secondary flows in the 
planes perpendicular to the ampoule axes. The flow visual- 
ization based on our numerical simulation shows a complex 
flow structure with multiple flow cells in the xy planes which 
vary along the ampoule axes. The typical pattern of the sec- 
ondary flow consists of two counterrotating flow cells on each 
side of the mid-center vertical (symmetric) plane (MCP). The 
upper flow cell flows down along the vertical rear (and front) 
wall and up at the MCP. The lower flow cell brings the fluid 
up along the vertical ampoule walls and down at the MCP. 

In the no-growth case, the 2-D model provides a good 
approximation for the interface shape and the primary flow 
pattern on the MCP. However, for non-zero growth rate, the 
2-D model is not adequate and 3-D modeling must be used. 
Comparisons between the two sets of solutions obtained from 
the round-comer and the simplified sharp-comer meshes also 
suggest that the shape of the ampoule comer does not have a 
significant effect on the numerical solution. 

Our experimental examination is continuing with 3-D 
fluid flow velocity measurements in the SCN liquid near the 
interface. The valuable experience gained from our com- 
parisons between experiments and simulations, and the good 
agreement achieved thus far, enables us to more confidently 
apply the code to other modelling cases (Yao et al„ 1993, Yao 
and de Groh, 1994). 
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